/*
 * Copyright 2010-2012 Susanta Tewari. <freecode4susant@users.sourceforge.net>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

package genomemap.stat;

import Jama.Matrix;

import genomemap.data.GDataSource;
import genomemap.worker.GLikeliHood;
import genomemap.worker.helper.StdErrHelper;

import javautil.io.LogUtil;
import javautil.lang.MathUtil;

/**
 * Class description
 *
 * @version        Enter version here..., 12/11/23
 * @author         Susanta Tewari
 */
public class StandardError {

    /** Field description */
    private static final double STD_CONVG_PRECISION = 1.0E-8;

    /**
     * @todo throw exception for possible errors
     * @param gHood
     * @param verbose
     * @return
     * @throws java.lang.Exception
     */
    public static double[] getStdErrs(GLikeliHood gHood, boolean verbose) throws Exception {

        /**
         * Compute the Information matrix for vector C
         * at cProbEst. This requires in turn computing
         * single and double derivatives at cProbEst and the computation
         * of the jacobian makes use of the iteration histroy.
         */
        GDataSource gData     = gHood.getGData();
        int[][] genotype      = gData.getArrayOfDistinctGenotypes(gData.getdata());
        int distinctGenotypes = genotype[0].length;
        int loci              = gData.getLoci();
        int[] genotypeFreq    = gData.getArrayOfCellFreq(gData.getdata());
        double[][] info = StdErrHelper.getInfoDerivs(genotype, genotypeFreq, gHood.getCProbEst(),
                              distinctGenotypes);

        if (verbose) {
            LogUtil.print(info);
        }

        Matrix delV    = null;
        Matrix inverse = null;

        /**
         * Compute the Jacobian
         */
        try {

            inverse = new Matrix(info).inverse();

            Matrix dm2Mat = getDm2Mat(gHood.getCProbEst(), gData, gHood.getIterationHistory(),
                                      verbose);

            delV = inverse.times(dm2Mat).times((Matrix.identity(loci - 1,
                    loci - 1).minus(dm2Mat)).inverse());

        } catch (Exception exception) {
            throw new Exception("Computation of Standard Error failed:");
        }

        Matrix infoMatrix = inverse.plus(delV);

        /**
         * Compute the standard errors of C using the information matrix
         */
        double[] stdErrs = new double[infoMatrix.getRowDimension()];

        for (int index = 0; index < infoMatrix.getRowDimension(); index++) {


            // System.out.println (infoMatrix.get(index, index));
            stdErrs[index] = Math.sqrt(infoMatrix.get(index, index));
        }

        return stdErrs;
    }

    /**
     * One co-ordinate of cProbFinal is tweaked and one EM iteration
     * is run to get the new cProbEst.
     *
     * @param cProbFinal description
     * @param gData description
     * @param cProbi description
     * @param index description
     * @param verbose description
     *
     * @return description
     */
    private static double[] RunEM(double[] cProbFinal, GDataSource gData, double cProbi, int index,
                                  boolean verbose) {


        // cProbFinal is tweaked
        double[] cProbOldTweaked = cProbFinal.clone();

        cProbOldTweaked[index] = cProbi;


        // 1 EM iteration is run on the tweaked parameter vector
        GLikeliHood gHood = new GLikeliHood();

        gHood.setCProbInitial(cProbOldTweaked);
        gHood.setVerbose(verbose);
        gHood.setIterationLimit(1);
        gHood.setConvgEpsilon(STD_CONVG_PRECISION);
        gHood.setGData(gData);
        gHood.runGHood();

        return gHood.getCProbEst();
    }

    /**
     * Method description
     *
     * @param cProbFinal description
     * @param gData description
     * @param iterates description
     * @param verbose description
     *
     * @return description
     */
    private static Matrix getDm2Mat(double[] cProbFinal, GDataSource gData, double[][] iterates,
                                    boolean verbose) {

        int loci       = cProbFinal.length + 1;
        int[] semCode  = new int[loci - 1];
        int[] semCode2 = new int[loci - 1];
        int countZero  = 0,
            countOne   = 0,
            countZero1 = 0,
            countOne1  = 0;
        int[] dmIndex  = new int[loci - 1];

        for (int index = 0; index < loci - 1; index++) {

            if (cProbFinal[index] == 0) {

                semCode[index]  = 0;
                semCode2[index] = 0;

                countZero++;

            } else {

                semCode[index]  = 1;
                semCode2[index] = 1;
                dmIndex[index]  = countOne;    // observe that dmIndex containes indices for dm rows.

                countOne++;
            }
        }

        int[] semZeroIndex = new int[countZero];
        int[] semOneIndex  = new int[countOne];

        for (int index = 0; index < loci - 1; index++) {

            if (cProbFinal[index] == 0) {

                semZeroIndex[countZero1] = index;

                countZero1++;

            } else {

                semOneIndex[countOne1] = index;

                countOne1++;
            }
        }

        double[][] dm1 = new double[loci - 1][loci - 1];
        double[][] dm2 = new double[loci - 1][loci - 1];


        // run the SEM iterations
        for (int semIteration = 0; semIteration < iterates.length; semIteration++) {

            for (int index1 = 0; index1 < loci - 1; index1++) {

                if ((semCode[index1] == 1) && (semCode2[index1] == 1)) {

                    double[] semValue = RunEM(cProbFinal, gData, iterates[semIteration][index1],
                                              index1, verbose);

                    for (int index2 = 0; index2 < loci - 1; index2++) {

                        if (semCode[index2] == 1) {

                            if (semIteration == 0) {

                                dm1[dmIndex[index1]][dmIndex[index2]] = (semValue[index2]
                                        - cProbFinal[index2]) / (iterates[semIteration][index1]
                                                                 - cProbFinal[index1]);

                            } else {

                                dm2[dmIndex[index1]][dmIndex[index2]] = (semValue[index2]
                                        - cProbFinal[index2]) / (iterates[semIteration][index1]
                                                                 - cProbFinal[index1]);
                            }
                        }
                    }


                    // check for convergence : if acheived set the semCode to 0 for that row
                    if (semIteration > 0) {


                        // i.e., the first iteration has passed
                        semCode2[index1] = 0;

                        for (int index2 = 0; index2 < loci - 1; index2++) {

                            if (semCode[index2] == 1) {

                                if (Math.abs(dm1[dmIndex[index1]][dmIndex[index2]]
                                             - dm2[dmIndex[index1]][dmIndex[index2]]) > 0.01) {
                                    semCode2[index1] = 1;
                                }
                            }
                        }
                    }
                }
            }

            if (verbose) {
                System.out.println("SEM iteration record at iteration = " + semIteration);
            }

            if (MathUtil.Sum(semCode2) == 0) {

                if (verbose) {
                    System.out.println("SEM has converged at iteration = " + semIteration);
                }

                break;    // all rows have converged

            } else {
                dm1 = dm2.clone();    // old and new are swapped
            }
        }

        return new Matrix(dm2);
    }
}    // ends the class

